/*************************************************************************
 *                                                                       *
 * This file is part of the DogosFrame distribution.  For full copyright *
 * information, see COPYRIGHT and COPYING.LESSER.                        *
 *                                                                       *
 * Author:           sqbang                                              *
 * Email:            sqbang@ustc.edu                                     *
 *                                                                       *
 *************************************************************************/

#ifndef DOGOS_include_COMMON_math_TypeVector_h
#define DOGOS_include_COMMON_math_TypeVector_h

//COMMON headers
#include "COMMON/base/Common.h"
#include "COMMON/tbox/CompareTypes.h"

//C++ headers
#include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
#include <cmath>
#include <iostream> //std::cout
#include <complex>

namespace DOGOS
{

// Forward declarations
template <typename T> class VectorValue;

/**
 * This class defines a vector in \p 3 dimensional space of type T.
 * T may either be double or Complex.  The default constructor for
 * this class is protected, suggesting that you should not instantiate
 * one of these directly.  Instead use one of the derived types: \p Point
 * for a real-valued point in 3-space, or \p SpaceVector for a real
 * or complex-valued vector in 3-space.
 *
 * \author Benjamin S. Kirk
 * \date 2003
 */
template <typename T>
class TypeVector
{
    template <typename T2>
    friend class TypeVector;

protected:

    /**
     * Empty constructor.  Gives the vector 0 in \p 3
     * dimensions.
     */
    TypeVector  ();


    /**
     * Constructor-from-T.  By default sets higher dimensional
     * entries to 0.
     */
    TypeVector (const T x,
                const T y=0,
                const T z=0);


    /**
     * Constructor-from-scalars.  By default sets higher dimensional
     * entries to 0.
     */
    template <typename Scalar1, typename Scalar2, typename Scalar3>
    TypeVector (typename
                boostcopy::enable_if_c<ScalarTraits<Scalar1>::value,
                const Scalar1>::type x,
                typename
                boostcopy::enable_if_c<ScalarTraits<Scalar2>::value,
                const Scalar2>::type y=0,
                typename
                boostcopy::enable_if_c<ScalarTraits<Scalar3>::value,
                const Scalar3>::type z=0);


    /**
     * Constructor-from-scalar.  Sets higher dimensional entries to 0.
     * Necessary because for some reason the constructor-from-scalars
     * alone is insufficient to let the compiler figure out
     * TypeVector<Complex> v = 0;
     */
    template <typename Scalar>
    TypeVector (const Scalar x,
                typename
                boostcopy::enable_if_c<ScalarTraits<Scalar>::value,
                const Scalar>::type * sfinae = DOGOS_NULL_PTR);

public:

    /**
     * Helper typedef for C++98 generic programming
     */
    typedef T value_type;

    /**
     * Copy-constructor.
     */
    template <typename T2>
    TypeVector (
        const TypeVector<T2> & rhs);

    /**
     * Destructor.
     */
    ~TypeVector ();

    /**
     * Assign to a vector without creating a temporary.
     */
    template <typename T2>
    void
    assign (
        const TypeVector<T2>& rhs);

    /**
     * Assignment-from-scalar operator.  Used only to zero out vectors.
     */
    template <typename Scalar>
    typename boostcopy::enable_if_c<
    ScalarTraits<Scalar>::value,
                 TypeVector &>::type
                 operator = (const Scalar & p)
    { TBOX_ASSERT_EQ (p, Scalar(0)); this->zero(); return *this; }

    /**
     * Return the \f$ i^{th} \f$ element of the vector.
     */
    const T & operator () (const unsigned int i) const;

    const T & slice (const unsigned int i) const { return (*this)(i); }

    /**
     * Return a writeable reference to the \f$ i^{th} \f$ element of the vector.
     */
    T & operator () (const unsigned int i);

    T & slice (const unsigned int i) { return (*this)(i); }

    /**
     * Add two vectors.
     */
    template <typename T2>
    TypeVector<typename CompareTypes<T, T2>::supertype>
    operator + (const TypeVector<T2> &) const;

    /**
     * Add to this vector.
     */
    template <typename T2>
    const TypeVector<T> & operator += (const TypeVector<T2> &);

    /**
     * Add to this vector without creating a temporary.
     */
    template <typename T2>
    void add (const TypeVector<T2> &);

    /**
     * Add a scaled value to this vector without
     * creating a temporary.
     */
    template <typename T2>
    void add_scaled (const TypeVector<T2> &, const T);

    /**
     * Subtract two vectors.
     */
    template <typename T2>
    TypeVector<typename CompareTypes<T, T2>::supertype>
    operator - (const TypeVector<T2> &) const;

    /**
     * Subtract from this vector.
     */
    template <typename T2>
    const TypeVector<T> & operator -= (const TypeVector<T2> &);

    /**
     * Subtract from this vector without creating a temporary.
     */
    template <typename T2>
    void subtract (const TypeVector<T2> &);

    /**
     * Subtract a scaled value from this vector without
     * creating a temporary.
     */
    template <typename T2>
    void subtract_scaled (const TypeVector<T2> &, const T);

    /**
     * Return the opposite of a vector
     */
    TypeVector<T> operator - () const;

    /**
     * Multiply a vector by a number, i.e. scale.
     */
    template <typename Scalar>
    typename boostcopy::enable_if_c<
    ScalarTraits<Scalar>::value,
                 TypeVector<typename CompareTypes<T, Scalar>::supertype> >::type
                 operator * (const Scalar) const;

    /**
     * Multiply this vector by a number, i.e. scale.
     */
    const TypeVector<T> & operator *= (const T);

    /**
     * Divide a vector by a number, i.e. scale.
     */
    template <typename Scalar>
    typename boostcopy::enable_if_c<
    ScalarTraits<Scalar>::value,
                 TypeVector<typename CompareTypes<T, Scalar>::supertype> >::type
                 operator / (const Scalar) const;

    /**
     * Divide this vector by a number, i.e. scale.
     */
    const TypeVector<T> & operator /= (const T);

    /**
     * Multiply 2 vectors together, i.e. dot-product.
     * The vectors may be of different types.
     */
    template <typename T2>
    typename CompareTypes<T, T2>::supertype
    operator * (const TypeVector<T2> &) const;

    /**
     * Multiply 2 vectors together, i.e. dot-product.
     * The vectors may be of different types.
     */
    template <typename T2>
    typename CompareTypes<T, T2>::supertype
    contract (const TypeVector<T2> &) const;

    /**
     * Cross 2 vectors together, i.e. cross-product.
     */
    template <typename T2>
    TypeVector<typename CompareTypes<T, T2>::supertype>
    cross(const TypeVector<T2> &) const;

    /**
     * Think of a vector as a \p dim dimensional vector.  This
     * will return a unit vector aligned in that direction.
     */
    TypeVector<T> unit() const;

    /**
     * Returns the magnitude of the vector, i.e. the square-root of the
     * sum of the elements squared.  This function is deprecated, instead
     * call norm().
     */
    double size() const;

    /**
     * Returns the magnitude of the vector, i.e. the square-root of the
     * sum of the elements squared.
     */
    double norm() const;

    /**
     * Returns the magnitude of the vector squared, i.e. the sum of the
     * element magnitudes squared.
     */
    double norm_sq() const;

    /**
     * Zero the vector in any dimension.
     */
    void zero();

    /**
     * @returns \p true iff two vectors occupy approximately the same
     * physical location in space, to within a relative tolerance of \p tol.
     */
    bool relative_fuzzy_equals(const TypeVector<T> & rhs,
                               double tol = TBOX_TOL) const;

    /**
     * @returns \p true iff two vectors occupy approximately the same
     * physical location in space, to within an absolute tolerance of \p tol.
     */
    bool absolute_fuzzy_equals(const TypeVector<T> & rhs,
                               double tol = TBOX_TOL) const;

    /**
     * @returns \p true iff this(i)==rhs(i) for each component of the
     * vector. For floating point types T, the function
     * absolute_fuzzy_equals may be a more appropriate choice.
     */
    bool operator == (const TypeVector<T> & rhs) const;

    /**
     * @returns \p true iff two vectors do not occupy approximately the same
     * physical location in space.
     */
    bool operator != (const TypeVector<T> & rhs) const;

    /**
     * @returns \p true if this vector is "less"
     * than another.  Useful for sorting.
     * Also used for choosing some arbitrary basis function
     * orientations
     */
    bool operator < (const TypeVector<T> & rhs) const;

    /**
     * @returns \p true if this vector is "less"
     * than or equal to another.  Useful for sorting.
     * Also used for choosing some arbitrary constraint
     * equation directions
     */
    bool operator <= (const TypeVector<T> & rhs) const;

    /**
     * @returns \p true if this vector is "greater"
     * than another.  Useful for sorting.
     * Also used for choosing some arbitrary basis function
     * orientations
     */
    bool operator > (const TypeVector<T> & rhs) const;

    /**
     * @returns \p true if this vector is "greater"
     * than or equal to another.  Useful for sorting.
     * Also used for choosing some arbitrary constraint
     * equation directions
     */
    bool operator >= (const TypeVector<T> & rhs) const;

    /**
     * Formatted print, by default to \p libMesh::out.
     */
    void print(std::ostream & os = std::cout) const;

    /**
     * Formatted print as above but allows you to do
     * Point p(1,2,3);
     * std::cout << p << std::endl;
     */
    friend std::ostream & operator << (std::ostream & os, const TypeVector<T> & t)
    {
        t.print(os);
        return os;
    }

    /**
     * Unformatted print to the stream \p out.  Simply prints the elements
     * of the vector separated by spaces.  Optionally prints a newline,
     * which it does by default.
     */
    void write_unformatted (std::ostream & out, const bool newline = true) const;

protected:

    /**
     * The coordinates of the \p TypeVector
     */
    T _coords[3];
};





//------------------------------------------------------
// Inline functions

template <typename T>
inline
TypeVector<T>::TypeVector ()
{
    _coords[0] = 0;
    _coords[1] = 0;
    _coords[2] = 0;
}



template <typename T>
inline
TypeVector<T>::TypeVector (const T x,
                           const T y,
                           const T z)
{
    _coords[0] = x;
    _coords[1] = y;
    _coords[2] = z;
}


template <typename T>
template <typename Scalar1, typename Scalar2, typename Scalar3>
inline
TypeVector<T>::TypeVector (typename
                           boostcopy::enable_if_c<ScalarTraits<Scalar1>::value,
                           const Scalar1>::type x,
                           typename
                           boostcopy::enable_if_c<ScalarTraits<Scalar2>::value,
                           const Scalar2>::type y,
                           typename
                           boostcopy::enable_if_c<ScalarTraits<Scalar3>::value,
                           const Scalar3>::type z)
{
    _coords[0] = x;
    _coords[1] = y;
    _coords[2] = z;
}



template <typename T>
template <typename Scalar>
inline
TypeVector<T>::TypeVector (const Scalar x,
                           typename
                           boostcopy::enable_if_c<ScalarTraits<Scalar>::value,
                           const Scalar>::type * /*sfinae*/)
{
    _coords[0] = x;
    _coords[1] = 0;
    _coords[2] = 0;
}



template <typename T>
template <typename T2>
inline
TypeVector<T>::TypeVector (const TypeVector<T2> & p)
{
    // copy the nodes from vector p to me
    for (unsigned int i=0; i<3; i++)
        _coords[i] = p._coords[i];
}



template <typename T>
inline
TypeVector<T>::~TypeVector ()
{
}



template <typename T>
template <typename T2>
inline
void TypeVector<T>::assign (const TypeVector<T2> & p)
{
    for (unsigned int i=0; i<3; i++)
        _coords[i] = p._coords[i];
}



template <typename T>
inline
const T & TypeVector<T>::operator () (const unsigned int i) const
{
    TBOX_ASSERT_LT (i, 3);

    return _coords[i];
}



template <typename T>
inline
T & TypeVector<T>::operator () (const unsigned int i)
{
    TBOX_ASSERT_LT (i, 3);

    return _coords[i];
}



template <typename T>
template <typename T2>
inline
TypeVector<typename CompareTypes<T, T2>::supertype>
TypeVector<T>::operator + (const TypeVector<T2> & p) const
{
    typedef typename CompareTypes<T, T2>::supertype TS;

    return TypeVector<TS> (_coords[0] + p._coords[0],
                           _coords[1] + p._coords[1],
                           _coords[2] + p._coords[2]);
}



template <typename T>
template <typename T2>
inline
const TypeVector<T> & TypeVector<T>::operator += (const TypeVector<T2> & p)
{
    this->add (p);

    return *this;
}



template <typename T>
template <typename T2>
inline
void TypeVector<T>::add (const TypeVector<T2> & p)
{
    _coords[0] += p._coords[0];
    _coords[1] += p._coords[1];
    _coords[2] += p._coords[2];
}



template <typename T>
template <typename T2>
inline
void TypeVector<T>::add_scaled (const TypeVector<T2> & p, const T factor)
{
    _coords[0] += factor*p(0);
    _coords[1] += factor*p(1);
    _coords[2] += factor*p(2);
}



template <typename T>
template <typename T2>
inline
TypeVector<typename CompareTypes<T, T2>::supertype>
TypeVector<T>::operator - (const TypeVector<T2> & p) const
{
    typedef typename CompareTypes<T, T2>::supertype TS;

    return TypeVector<TS>(_coords[0] - p._coords[0],
                          _coords[1] - p._coords[1],
                          _coords[2] - p._coords[2]);
}



template <typename T>
template <typename T2>
inline
const TypeVector<T> & TypeVector<T>::operator -= (const TypeVector<T2> & p)
{
    this->subtract (p);

    return *this;
}



template <typename T>
template <typename T2>
inline
void TypeVector<T>::subtract (const TypeVector<T2> & p)
{
    for (unsigned int i=0; i<3; i++)
        _coords[i] -= p._coords[i];
}



template <typename T>
template <typename T2>
inline
void TypeVector<T>::subtract_scaled (const TypeVector<T2> & p, const T factor)
{
    for (unsigned int i=0; i<3; i++)
        _coords[i] -= factor*p(i);
}



template <typename T>
inline
TypeVector<T> TypeVector<T>::operator - () const
{
    return TypeVector(-_coords[0],
                      -_coords[1],
                      -_coords[2]);
}



template <typename T>
template <typename Scalar>
inline
typename boostcopy::enable_if_c<
ScalarTraits<Scalar>::value,
             TypeVector<typename CompareTypes<T, Scalar>::supertype> >::type
             TypeVector<T>::operator * (const Scalar factor) const
{
    typedef typename CompareTypes<T, Scalar>::supertype SuperType;

    return TypeVector<SuperType>(_coords[0]*factor,
                                 _coords[1]*factor,
                                 _coords[2]*factor);
}



template <typename T, typename Scalar>
inline
typename boostcopy::enable_if_c<
ScalarTraits<Scalar>::value,
             TypeVector<typename CompareTypes<T, Scalar>::supertype> >::type
             operator * (const Scalar factor,
                         const TypeVector<T> & v)
{
    return v * factor;
}



template <typename T>
inline
const TypeVector<T> & TypeVector<T>::operator *= (const T factor)
{
    _coords[0] *= factor;
    _coords[1] *= factor;
    _coords[2] *= factor;

    return *this;
}



template <typename T>
template <typename Scalar>
inline
typename boostcopy::enable_if_c<
ScalarTraits<Scalar>::value,
             TypeVector<typename CompareTypes<T, Scalar>::supertype> >::type
             TypeVector<T>::operator / (const Scalar factor) const
{
    TBOX_ASSERT_NE (factor, static_cast<T>(0.));

    typedef typename CompareTypes<T, Scalar>::supertype TS;

    return TypeVector<TS>(_coords[0]/factor,
                          _coords[1]/factor,
                          _coords[2]/factor);
}




template <typename T>
inline
const TypeVector<T> &
TypeVector<T>::operator /= (const T factor)
{
    TBOX_ASSERT_NE (factor, static_cast<T>(0.));

    for (unsigned int i=0; i<3; i++)
        _coords[i] /= factor;

    return *this;
}




template <typename T>
template <typename T2>
inline
typename CompareTypes<T, T2>::supertype
TypeVector<T>::operator * (const TypeVector<T2> & p) const
{
    return (_coords[0]*p(0) +
            _coords[1]*p(1) +
            _coords[2]*p(2));
}

template <typename T>
template <typename T2>
inline
typename CompareTypes<T, T2>::supertype
TypeVector<T>::contract(const TypeVector<T2> & p) const
{
    return (*this)*(p);
}



template <typename T>
template <typename T2>
TypeVector<typename CompareTypes<T, T2>::supertype>
TypeVector<T>::cross(const TypeVector<T2> & p) const
{
    typedef typename CompareTypes<T, T2>::supertype TS;

    // |     i          j          k    |
    // |(*this)(0) (*this)(1) (*this)(2)|
    // |   p(0)       p(1)       p(2)   |

    return TypeVector<TS>( _coords[1]*p._coords[2] - _coords[2]*p._coords[1],
                           -_coords[0]*p._coords[2] + _coords[2]*p._coords[0],
                           _coords[0]*p._coords[1] - _coords[1]*p._coords[0]);
}



template <typename T>
inline
double TypeVector<T>::size() const
{
    TBOX_DEPRECATED();
    return this->norm();
}



template <typename T>
inline
double TypeVector<T>::norm() const
{
    return std::sqrt(this->norm_sq());
}


template <typename T>
inline
double TypeVector<T>::norm_sq() const
{
    return (_coords[0]*_coords[0] +
            _coords[1]*_coords[1] +
            _coords[2]*_coords[2]);
}


template <typename T>
inline
void TypeVector<T>::zero()
{
    for (unsigned int i=0; i<3; i++)
        _coords[i] = 0.;
}


template <typename T>
inline
bool TypeVector<T>::absolute_fuzzy_equals(const TypeVector<T> & rhs, double tol) const
{
    return (std::abs(_coords[0] - rhs._coords[0]) +
            std::abs(_coords[1] - rhs._coords[1]) +
            std::abs(_coords[2] - rhs._coords[2])
            <= tol);
}



template <typename T>
inline
bool TypeVector<T>::relative_fuzzy_equals(const TypeVector<T> & rhs, double tol) const
{
    return this->absolute_fuzzy_equals(rhs, tol *
                                       (std::abs(_coords[0]) + std::abs(rhs._coords[0]) +
                                        std::abs(_coords[1]) + std::abs(rhs._coords[1]) +
                                        std::abs(_coords[2]) + std::abs(rhs._coords[2])));
}



template <typename T>
inline
bool TypeVector<T>::operator == (const TypeVector<T> & rhs) const
{
    return (_coords[0] == rhs._coords[0] &&
            _coords[1] == rhs._coords[1] &&
            _coords[2] == rhs._coords[2]);
}



template <typename T>
inline
bool TypeVector<T>::operator != (const TypeVector<T> & rhs) const
{
    return (!(*this == rhs));
}


//------------------------------------------------------
// Non-member functions on TypeVectors

// Compute a * (b.cross(c)) without creating a temporary
// for the cross product.  Equivalent to the determinant
// of the 3x3 tensor:
// [a0, a1, a2]
// [b0, b1, b2]
// [c0, c1, c2]
template <typename T>
inline
T triple_product(const TypeVector<T> & a,
                 const TypeVector<T> & b,
                 const TypeVector<T> & c)
{
    return
        a(0)*(b(1)*c(2) - b(2)*c(1)) -
        a(1)*(b(0)*c(2) - b(2)*c(0)) +
        a(2)*(b(0)*c(1) - b(1)*c(0));
}



/**
 * Compute |b x c|^2 without creating the extra temporary produced by
 * calling b.cross(c).norm_sq().
 */
template <typename T>
inline
T cross_norm_sq(const TypeVector<T> & b,
                const TypeVector<T> & c)
{
    T x = b(1)*c(2) - b(2)*c(1),
      y = b(0)*c(2) - b(2)*c(0),
      z = b(0)*c(1) - b(1)*c(0);

    return x*x + y*y + z*z;
}



/**
 * Calls cross_norm_sq() and takes the square root of the result.
 */
template <typename T>
inline
T cross_norm(const TypeVector<T> & b,
             const TypeVector<T> & c)
{
    return std::sqrt(cross_norm_sq(b,c));
}

} // namespace DOGOS

#include "COMMON/math/TypeVector.inl"

#endif
